Consistency of post-Newtonian waveforms with numerical relativity 
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General relativity predicts the gravitational wave signatures of coalescing binary black holes. Ex- 
plicit waveform predictions for such systems, required for optimal analysis of observational data, 
have so far been achieved using the post-Newtonian (PN) approximation. The quality of this treat- 
ment is unclear, however, for the important late-inspiral portion. We derive late-inspiral waveforms 
via a complementary approach, direct numerical simulation of Einstein's equations. We compare 
waveform phasing from simulations of the last ~ 14 cycles of gravitational radiation from equal- 
mass, nonspinning black holes with the corresponding 2.5PN, 3PN, and 3.5PN orbital phasing. We 
find phasing agreement consistent with internal error estimates based on either approach, suggesting 
that PN waveforms for this system are effective until the last orbit prior to final merger. 

PACS numbers: 04.25.Dm, 04.25.Nx 04.30.-w 04.30.Db, 95.30.Sf, 97.60.Lf 



Compact astrophysical binaries spiral together due to 
the emission of gravitational radiation. Calculating the 
dynamics of these systems, and the corresponding gravi- 
tational waveforms, has been a central problem in general 
relativity for several decades. With first-generation inter- 
ferometers such as LIGO, VIRGO, and GEO600 now op- 
erating, and development moving forward on the space- 
based LISA mission, accurate and reliable waveforms are 
urgently needed for gravitational-wave data analysis. 

Post-Newtonian (PN) methods, based on expansions 
in the parameter e ~ v/c, have been the major ana- 
lytic tool used to calculate the system dynamics and 
waveforms during the early part of the inspiral, when 
the binary components are relatively widely separated 
and thus have a small orbital frequency Currently, 
gravitational-wave data analysis for binary inspiral relies 
on waveforms derived from PN methods [2J • The current 
predicted orbital phase is available up to 0(e 7 ), which 
is referred to as 3.5PN order. However the convergence 
properties of the PN sequence are not well understood, 
and it is not yet clear how well PN predictions work late 
in the inspiral when frequencies arc high. 

Numerical relativity, in which the full set of Einstein's 
equations is solved on a computer, is needed to handle 
the final stages of the binary evolution, when the com- 
ponents inspiral rapidly and merge. Recently, there has 
been dramatic progress in the use of numerical relativity 
to simulate the final insp iral and merger of black holes 
0, II S 0, 0, S, H M, Ef • These breakthroughs have al- 
lowed numerical simulations with increasingly wider ini- 
tial separations, producing longer wavetrains. Linking 
such simulations with the PN calculations and comparing 
their waveforms in the late inspiral regime is a pressing 
concern of gravitational- wave data analysis. While quali- 
tative comparisons have suggested that the two methods 
agree fairly well until shortly before the merger 0, EH , 



our group has found that even inaccurate waveforms may 
approximately agree over several cycles. Quantifying the 
level of agreement unambiguously requires long, accurate 
simulations with very low eccentricity, and careful anal- 
ysis. 

We have carried out suitable numerical simulations of a 
merging equal-mass, nonspinning black-hole binary. The 
black holes start on nearly circular orbiting trajectories 
~ 1200A/ before merger, where Af refers to the mass that 
the system would have had when the black holes were still 
far apart, before radiative losses were significant. M is 
related to time by M = 5x lfT~ 6 s(M /A/©). In this Letter, 
we quantitatively compare crucial phasing information in 
our numerical simulation waveforms with phasing in PN 
waveforms, finding striking agreement. 

The numerical simulations were performed using the 
moving puncture method 0, 0, E2] • We use fourth-order 
Runge-Kutta time integration, fourth-order-accurate fi- 
nite spatial differencing, and second-order-accurate ini- 
tial data. Adaptive mesh refinement is used to resolve 
both the dynamics near the black holes and the propaga- 
tion of the gravitational waves 0|. The wave extraction 
sphere was of radius 60M and the cubical outer boundary 
was of half-width 1536M; more details about this simu- 
lation can be found in [1 31 ] - We performed physically 
equivalent runs at three different maximum resolutions: 
low (3A//64), medium (3A//80), and high (Af/32). We 
find fourth-order convergence of the Hamiltonian con- 
straint, and better than second-order convergence of the 
momentum constraints during the runs. 

The simulations begin at an angular gravitational-wave 
frequency ui ~ 0.051A/ -1 . The frequency then sweeps 
upward through roughly an order of magnitude while the 
black holes undergo ~ 7 orbits, producing ~ 14 gravita- 
tional wave cycles before merger. For such long-lasting 
simulations, the primary consideration in providing a re- 
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FIG. 1: Gravitational strain waveforms from the merger of 
equal-mass Schwarzschild black holes. The solid curve is the 
waveform from the high resolution numerical simulation, and 
the dashed curve is a PN waveform with 3.5PN order phasing 
[3, EH and 2.5PN order amplitude accuracy plj. Time t — 
is the moment of peak radiation amplitude in the simulation. 



alistic initial data model is to set up the orbiting black 
holes with minimal eccentricity, as gravitating binary sys- 
tems of comparable-mass objects are expected to circu- 
larize rapidly through the emission of gravitational radia- 
tion. We have selected an initial black hole configuration 
with the relatively low eccentricity of less than one per- 
cent, as measured below. 

Fig.[T]shows the gravitational wave strain generated by 
our highest-resolution numerical run and that predicted 
by the PN approximation with 3.5PN phasing [ij], [HI 
and 2.5PN (beyond leading order) amplitude accuracy 
[lij ]. The waves are based on the dominant I = 2, m = 2 
spin-weighted spherical harmonic of the radiation, and 
represent an observation made on the system's equato- 
rial plane, where only one polarization component con- 
tributes to the measured strain. The initial phase and 
time of the waves have been adjusted so that the fre- 
quency and phase for each waveform agree at a point, 
t = — 1000M, that is early in the simulation, but af- 
ter transient effects from the initial data have subsided. 
We will quantify the phase agreement below using the 
frequency domain, so that the time shifting, done for il- 
lustrative purposes in Fig. [1] will have no impact on the 
subsequent analysis. 

To conduct comparisons with PN calculations, wc need 
to extract an instantaneous gauge-invariant polarization 
phase 4> and angular frequency to from our simulations. 
These are derived from the gravitational wave strain's 
first time-derivative, which is a robust quantity in the 
numerical data. This frequency corresponds to the sweep 
rate of the polarization angle of the circularly polarized 



gravitational wave that can be observed on the system's 
rotation axis. 

We define eccentricity as a deviation from an under- 
lying smooth, secular trend. We obtain a monotonic 
"secular" frequency-time relation by modeling the wave- 
form angular frequency u> as a fourth-order monotonic 
polynomial Lu c (t), plus an eccentric modulation of the 
form du)(t) = uj{t) - uj c (t) = Asin($(t)), where &(t) is 
a quadratic function of time. Fitting this equation to 
our data yields A = 8(±1) x 10~ 4 M _1 . For Keplerian 
systems, conserved angular momentum is proportional to 
r 2 oj, so the eccentricity corresponds to half the fractional 
amplitude of the frequency modulation: e = A/(2lu). In 
our case the eccentricity starts near 0.008, decreasing by a 
factor of three by the time oj c M ~ 0.15. We will compare 
our simulation with non-eccentric PN calculations, with 
the expectation that small eccentricities have minimal ef- 
fect on the important underlying secular trend in the rate 
at which frequency sweeps up approaching merger. 

The phasing of the waveform is critical for gravita- 
tional wave observation. For data analysis, the optimal 
methods for both detection and parameter estimation 
rely on matched filtering, which employs a weighted in- 
ner product that can be expressed in Fourier space as 
< h,s>=fdf[h*(f)s(f) + h(f)s*(f)]/S n (f), where h is 
the template being used, s is the signal being analyzed, 
and S n is the one-sided power spectral density of the de- 
tector's noise 17[. A template that maximizes < h,s > 



will provide an optimal filter. Therefore, the most crucial 
factor is the relative phasing of the template and signal. 
The inner product will cease to accumulate in sweeping 
through frequency if the template and the signal evolve 
to be out of phase with each other by more than a half- 
cycle, decreasing the effectiveness of the procedure. 

Our key objective is to compare phasing between nu- 
merical and PN waveforms. We can make a stronger con- 
nection to the underlying physics while avoiding issues 
with time alignment by comparing phases as a function 
of polarization frequency, which corresponds to twice the 
orbital frequency in the PN case. For circular inspiral 
this frequency should grow monotonically in time, with 
the frequency lu c providing a physical reference of the 
"hardness" of the tightening binary. 

Circular inspiral phasing information is typically de- 
rived in PN theory by imposing an energy balance re- 
lation to deduce the rate at which oj c evolves from the 
radiation rate at a specified value of u c [H . Though not 
strictly derived in the PN context, this physically sen- 
sible condition currently allows the determination of the 
chirp rate uj c (lu c ) up to 3.5PN order From such a re- 
lation, information about phase and time are determined 
by integrating d<f)/doj c — u c /lu c and dt/duj c = 1/lu c . The 
phasing information can be represented by any one of sev- 
eral relations between phase, frequency and time. Vari- 
ous approaches take the PN-expanded representation of 
one of these relations as the PN "result" for waveform 
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FIG. 2: Gravitational wave phase, in radians, for numerical 
and PN waveforms. The solid curve is a Richardson extrap- 
olation of the numerical results. The solid curve agrees well 
with the phase obtained by numerically integrating the 3.5PN 
expansion of the chirp rate lo c (uJc)- Each successive PN order 
shown agrees better with the Richardson-extrapolated result, 
although this is not true of all the preceding terms in the PN 
sequence, since the sequence does not converge monotonically. 



phasing [l], [ll], [l8| ■ It has been demonstrated 11 1 that 
the PN expansion of lu c (lu c ), numerically integrated as 
needed, has the greatest utility for conducting compar- 
isons of phasing with numerical results during the late 
inspiral, and we adopt that convention. 

For the purpose of comparison with our numerical sim- 
ulations, we invert the monotonic function u) c (t) to obtain 
the phase as a function of frequency: 4>(lu c ) — (f>(t(u c )). 
Note that the effect of eccentricity is not removed from 
c6, though the "circularized" frequency oj c does provide 
the abscissa according to which phases are compared in 
the different treatments. 

Fig. [2] shows the wave phases as a function of uj c M; 
here each phase is adjusted by addition of a constant 
so that it vanishes at lu c M = 0.054, corresponding to 
the time 1000M before the radiation merger peak in the 
high-resolution numerical simulation. As demonstrated 
in Fig. [31 the sequence of numerical results converges at 
fourth order, allowing us to obtain a fifth-order-accurate 
result by Richardson extrapolation (thick solid line in 
Fig. E}. 

In Fig. [2] we compare the numerical simulation results 
for the phase with the numerically integrated PN expan- 
sion of the chirp rate at 2.5PN, 3PN, and 3.5PN order. 
The agreement of the extrapolated numerical result with 
the integrated PN chirp rate improves with each succes- 
sive order, with the 3.5PN result showing striking agree- 
ment up to about oj c M ~ 0.15. 

We look more quantitatively at this agreement in Fig. [3] 
by plotting the phase differences 5<p that accumulate be- 
tween different phase approximations, as a function of 
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FIG. 3: Gravitational wave phase error estimates. Differ- 
ences between phasing from the integrated 3.5PN chirp rate 
and Richardson extrapolation from the numerical simulations 
(solid curve) are small, and are consistent with internal error 
estimates for the numerical simulation results and the PN se- 
quence. Curves which involve numerical phases are smoothed 
to remove high frequency noise. 



frequency. The thick dash-dotted curve shows the phase 
differences between our medium- and high-resolution re- 
sults, while the thin dash-dotted curve shows the dif- 
ferences between our low- and medium-resolution re- 
sults, scaled so that for fourth-order convergence the 
curves should superpose. This is indeed observed to 
good approximation. A good estimate for the error of 
the phase in the high-resolution run is given by its dif- 
ference from the phase obtained by Richardson extrapo- 
lation; this comes out to ~ 93% of the med-high curve 
shown in the figure. Note that the cumulative errors in 
the numerically-generated waveforms accrue primarily at 
lower frequencies, scaling approximately as u>~ 5 ; the thin 
dashed curve shows a fit to the med-high curve with this 
scaling. The deviations of the med-high curve from this 
fit show the effect of eccentricity in our simulations. This 
accumulation of phase errors at lower frequencies makes 
sense generally since the simulations spend longer in that 
regime. 

Without monotonic convergence between the 2PN, 
2.5PN and 3PN at the frequencies considered here, it is 
difficult to estimate errors in the PN phase. Nonetheless 
we tentatively take the difference between the integrated 
3PN and 3.5PN chirp rates, shown by the dashed line in 
Fig. [3l as an upper bound on errors in a 3.5PN waveform. 
We also show (thin dotted line) the accumulated phase 
differences between the integrated 2.5PN and 3PN chirp 
rates to show the relative contribution of the preceding 
PN term. 

The trend in the slope of these error curves indicates 
the rate at which phase error accumulates, as indepen- 
dently estimated within each approach. The slope of 
the fit to the med-high error curve in Fig. [3] is ini- 
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tially higher than that of the 3PN-3.5PN error curve, 
but decreases steadily, matching the PN error-curve slope 
around uj c M ~ 0.08 and less thereafter. For clarity a 
piece of the 3PN-3.5PN curve has been translated up- 
ward to fit to the med-high curve in Fig. [3] This suggests 
that phasing errors for our high-resolution simulation ac- 
cumulate more quickly than 3.5 PN phasing errors for 
lu c M < 0.08 (t < -300M), with the numerical simu- 
lation phasing being more accurate than PN at higher 
frequencies. In both cases, the phase error accumulates 
to roughly two radians by u> c M ~ 0.15 as the black holes 
begin to plunge together. 

We now address the central objective of this Letter, 
a quantitative comparison of numerical and PN phasing 
results. We compare the Richardson-extrapolated phase, 
our best simulation estimate, with the integrated 3.5PN 
chirp rate result. The difference is shown in the solid 
curve of Fig. [3l Note that, overall, the phase differ- 
ences between the PN and numerical results (solid curve) 
accumulate less rapidly than the estimated errors in ci- 
ther the high resolution simulation or the 3.5PN results 
over much of the frequency range shown. This illustrates 
that the numerical and post-Newtonian predictions ap- 
pear to be converging to a common answer for waveform 
phase, in a frequency regime where the validity of the 
PN predictions could not previously be assessed. In the 
range 0.1 < lo c < 0.15 the slope of the trend in the solid 
curve surpasses that of the 3PN-3.5PN curve, with the 
precise crossover point obscured by the effect of eccen- 
tric oscillations in the numerical simulations; we have 
translated a portion of the 3PN-3.5PN curve down to 
the solid line for clarity. After this range, corresponding 
to —170 < t/M < —70, and occurring 1 to 3 wave cy- 
cles before the estimated end of the binary's last orbit, 
the phase difference between the best estimates of the 
numerical and PN approaches grows more significantly. 

Through the frequency range 0.054 < u> c M < 0.15 
the net phase difference, measured against frequency, 
amounts to less than one radian, a level of error that 
would be tolerable in many gravitational wave data anal- 
ysis applications. Superficially, phase differences in Fig.[T] 
may appear smaller than the differences quantified in 
Fig. [3l This is because the differences in chirp-rate are 
more directly evident in the frequency-based phase com- 
parisons. This more immediate connection to the merger 
dynamics, together with the avoidance of time-alignment 
issues, makes frequency-based phase comparisons a more 
reliable indicator of phasing differences. 

Our results provide a crucial cross-validation of PN 
waveforms from the late inspiral of binary merger, with 
results of new long-lasting numerical simulations. These 
simulations have sufficient accuracy to provide a mean- 
ingful comparison with PN waveforms over the last t ~ 
1000 M of the coalescence, specifically addressing a bi- 
nary system of equal-mass non-spinning black holes. We 
find phase agreement consistent with internal phase-error 



estimates conducted in each approach, indicating that 
phase accuracies within a few radians are now achievable 
for this part of the coalescence waveform. 

We emphasize, however, that there is still much impor- 
tant work to be done in improving and further assessing 
PN and numerical simulation waveforms. Certainly we 
have only addressed one case in a large parameter space 
of potential binaries, which will inevitably include sys- 
tems, such as rapidly precessing unequal-mass spinning 
binaries, that are harder to treat with present PN and nu- 
merical techniques. With either approach, even for our 
simple case, a non-negligible amount of phasing error ac- 
cumulates over the range studied, and more will have ac- 
cumulated at lower frequencies addressable through the 
PN approximation. We expect continuing developments 
in numerical simulations and the pursuit of higher-order 
PN treatments to be crucial for developing a refined un- 
derstanding of coalescence waveforms, which will be cru- 
cial in some data analysis applications for gravitational 
wave observations. 
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